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Abstract 

We consider the ordered and disordered dynamics for monolayers of rolling self-interacting par- 
ticles with an off-set center of mass and a non-isotropic inertia tensor. The rolling constraint is 
considered as a simplified model of a very strong, but rapidly decaying bond with the surface, 
preventing application of the standard tools of statistical mechanics. We show the existence and 
nonlinear stability of ordered lattice states, as well as disturbance propagation through and chaotic 
vibrations of these states. We also investigate the dynamics of disordered gas states and show that 
there is a surprising and robust linear connection between distributions of angular and linear ve- 
locity for both lattice and gas states, allowing to define the concept of temperature. 
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Introduction Molecular monolayers are playing an ever increasing role in technology as 
they allow manipulation of contact properties of materials in a precise and controlled man- 
ner. There have been considerable effort in analyzing the structure of molecular monolayers; 
in particular, self-assembled polymeric monolayers jTJ [2] and water monolayers with TIP5P 
model [3]. Previous theoretical work on molecular monolayers has concentrated on either 
static studies of monolayer structures [U [4] or molecular dynamics simulations where each 
molecule is moving under the forces and torques from surrounding molecules and the sub- 
strate [2, El El |6] . The direct molecular simulations employed in these papers were quite 
successful in explaining various properties of liquid water, e.g formation of contact angle [2] 
or anomalous properties for nano-confined water [5l|6]. In the latter papers, the confinement 
of water molecules to a very thin layer was achieved by physically constraining the bulk of 
water to a very narrow (nm) layer by the rigid surface on both sides. On the other hand, ex- 
perimental data demonstrate the presence of water monolayer on a Si surface under normal 
conditions, which is due to the strong bond between the water molecules and the substrate 
[HE]. Such a bond is difficult to account for in the atomistic molecular simulations of TIP5P 
type. On the other hand, one would like to keep the relative simplicity of atomistic models 
without introducing extra coupling bonds, the nature of which is not well understood. 
The purpose of this Letter is to suggest a simple yet physical way of introducing a strong 
coupling bond between the monoloayer molecules and the substrate as a rolling constrant. 
We assume that the bonds between the substrate and the molecule are very short ranged and 
thus act on the part of the molecule in close contact with the substrate atoms. Physically, 
that leads to the fact that while the molecule itself is moving, its point of contact with the 
substrate is stationary. This invokes the analogy with a classical problem of rolling body on 
the surface with perfect friction, appearing when the bond molecule-substrate is infinitely 
strong at contact point, but decays rapidly away from substrate. In reality, if that bond 
is large but finite, the motion will be a combination of sliding and rotating. However, the 
theory of sliding and rolling is not yet well developed [HI |9] , and so in this paper we will 
restrict ourselves to the simplest possible realization when the friction is perfect. Thus, we 
consider the dynamics of a monolayer consisting of rolling molecules that are self-interacting 
by the long-range interactions (Lennard- Jones and electrostatic), while the influence of the 
boundary is limited to restricting the motion to the perfect rolling dynamics. While the 
rolling dynamics may seem too idealized, it has actually been observed in the context of 
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functional nano-structures [ID]. More complex models of molecule-substrate interaction are 
possible, resulting in the gravity-like component in the force acting on the ball, but such 
interactions will not be considered here. 

The study of rolling motion of rigid bodies has a long history in the context of classical 
mechanics j9j HU H2]. However, the study of collective motion of rolling particles has not 
been undertaken. The non-holonomic rolling constraint is a major obstacle in the way of 
constructing statistical mechanics for the rolling particle systems | [T3l - FT6] . This Letter is 
devoted to defining ordered and disordered states of the rolling particle systems (akin to 
solid and gas/liquid), and computing statistical physics concepts for such systems. We also 
show the way rolling affects propagation of phonons through a rolling particle lattice, which 
opens the way to experimentally confirm the molecular rolling motion. 
Setup of the dynamics The rolling particles are simulated as identical spherical rigid 
bodies of radius r all having the same mass m and the moments of inertia tensors. The center 
of mass is assumed to be at a position different from the geometric center, as illustrated 
on Fig. [IJ The notation used in this paper is as follows: I\ is the unit vector pointing to 
the geometric center (GC), ft is the angular velocity in the ball's coordinate frame and is 
the vector pointing to the center of mass (CM). The equations of motion for an individual 



ball come from the well-known Chaplygin's equations j9j [IT]. Under cetrain symmetry 
conditions, Chaplygin's ball is completely integrable has three integrals of motion: one is 
total energy that is the easily understandable physical quantity. The other two constants 
of motion are Routh and Jellet integrals that are harder to explain in terms of elementary 
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FIG. 1: Schematic and definitions of an offset rolling ball dynamics. 
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physics [32]. The equations of motion for the i-th ball are as follows: 



d 
dt 



- mV x (fi* x s*) = x F* + T* 



(1) 




(2) 



Here, P are the eigenvalues of the tensor of inertia and F l and T l are the total force 
and torque acting on the zth particle. These forces and torques include the inter-particle 
interactions. The rolling constraint for velocity of the center of mass V*, written as V* = 
ft 1 x s l cannot be reduced to an equation between configuration variables only, and is thus 
non-holonomic [llj. The Hamiltonian approach for the ensemble of rolling ball cannot be 
defined, and thus it is impossible to construct a meaningful approach for this system by 
simply extending the classical statistical physics to this case. Note that neither linear nor 
angular momenta are conserved, neither for individual particles nor for the whole system. 
The total energy of the system, however, is conserved. 

The inter-molecular interactions causes the interchange of linear and angular momentum 
as well as energy between particles. In the simulation of rolling water molecules, the Lennard- 
Jones (LJ) potential Vlj and the dipole potential Vd are considered. The dipole is 
positioned at the center of mass on each spherical particle. The interactions we consider are 
defined as follows: 



To be concrete, in the computations of collective dynamics for rolling particles, we use 
the parameters and interactions relevant to water molecule monolayers by choosing the mass 
m = 2.991 • 1(T 23 g, moments of inertia (I u I 2 ,h) = (0.2076,0.1108,0.3184) • 10~ 39 g-cm 2 , 
radius r = 1 A, displacement of center of mass from the geometric center £ = 0.068 A, 
dipole moment 6.17 • 10 -30 (C • m), LJ radius a = 3.165 A and energy e = 0.650kJ/mol. 
These values correspond to the parameters of a water molecule [19j. For convenience, we 
choose the angular velocity scale Co = 10 13 rad/s. The rolling constraint then introduces the 
scaling of velocity to be v = ruj = lOcm/s. Given a different set of parameters, the details 
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of our computations will be different, but the methods and results outlined here hold for 
other molecules as well. For the set of parameters considered here, neither Routh nor Jellet 
integrals for each ball are conserved. 

Stationary states: a crystalline lattice The existence of stationary states for the system 
of rolling particles depends on the presence and orientation of the dipole moment. Suppose 
for now that the dipole moment is absent, and the only interaction between the molecules 
is Lennard- Jones. Suppose also for the moment that the particles are not moving, and are 
aligned so both geometric center and center of mass are along T, the unit vector pointing 
upward from the contact point. It is easy to see then that there is an equilibrium config- 
uration so that the centers of mass are arranged at a distance close to a, the equilibrium 
distance of LJ potential. For the case of two and three particles, there is an equilibrium 
configuration where the particles are arranged at exactly the distance a from each other. 

Let us suppose now that these particles are spun with angular velocity Hi = (0, 0, f^) 
that is pointing upwards. In the absence of nonlocal interaction this will be a neutrally 
stable state for each particle, although friction forces may destabilize such equilibria, similar 
to familiar phenomenon of the tippe top [20]. Because LJ forces are exactly at balance, 
independent of the values of f2 i3 this will be an equilibrium state. In the presence of the 
dipole moment, finite lattices cease to be equilibrium configurations. However, if the dipole 
moment points exactly along the line from the geometric center to the center of mass (which 
is the case for the particles considered here), an infinite regular lattice, possessing high 
symmetry (like a triangular or square lattice) will still be an equilibrium state. The dipole 
moments can either be arranged in the same directions or be alternating. However, the 
states with the dipole moments all pointing in one direction are unstable both linearly and 
nonlinear ly. Thus, in what follows we concentrate on the states with the alternating dipole 
moments as shown on Figure [2| 



The question about the stability of these states immediately arises. While a detailed 
investigation of the stability is a complex issue beyond the scope of this Letter, a simple 
physical explanation can be made showing that these states are linearly unstable. Suppose 




FIG. 2: Schematic representation of stationary states. 
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for simplicity that all rotation rates are the same, = Qq. A small perturbation in the 
spinning stationary state will result to a precession of each particle, which according to 
the rolling constraint will happen at the same rate as the rotation frequency. The distance 
between particles will also change with the same frequency, and through the LJ interaction 
([3]) both particles will experience parametric resonance. Similar argument applies to different 
rotation rates, and more particles in a lattice. Thus, all lattice states are linearly unstable. 
However, these states are nonlinearly stable at least for some of the configurations and 
initial conditions we have investigated, corresponding to low initial energies. In Fig. [3j left, 
we show the positions of the centers of 81 rolling balls situated in a rectangular lattice, over 
a long simulation time. While each particle remains close to its equilibrium position, the 
individual trajectories (blow-up on the same Figure) are chaotic and are strongly reminiscent 
of thermal vibrations in lattices. 



FIG. 3: Left: trajectory of the center of rolling particles in lattice model, with the blow-up showing 
the trajectory of an individual particle. The finite lattice is not an exact stationary state, so the 
vibrations are most apparent at the edges. Right: trajectories of a gas state in a circular container. 

Lattice dispersion relation In order to suggest a possible experimental verification of 
rolling motion, we suggest measuring the propagation of disturbances through the lattice. 
The derivation of dispersion relation for a lattice of spinning particles in general case is 
difficult, as it involves linearization about oscillating base states. Here, we present the 
analysis of the rolling and rocking lattice with all the particles being stationary in the base 
state, i.e. = 0. Assuming an infinite square lattice as described above, we arrive to the 
following dispersion relation for the propagation of disturbances of the form e - % ut+ik x x+ik v y ^ 
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(k x ,k y ) being the wave vector: 




(5) 
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where Q = /$/ (m(r + £) 2 ), K = cPVlj / dr 2 is the spring constant of the LJ potential and a 
is the periodicity of the square lattice . Note that ^ differs from the standard dispersion 
relation for a square lattice of springs only by the dimensionless coefficients incorporating 
the effects of rolling. For our values of parameters, Q ~ 0.1 — 0.2. Thus, the rolling constraint 
affects the speed of sound by about 10-20%, which should be a measurable difference. 
Disordered states: statistical analysis For large initial energies, the lattices become 
unstable, and in the absence of external boundaries the particles scatter to infinity, making 
the concepts of statistical physics meaningless. Thus, for large energies leading to gaseous 
states, we perform the simulations in a round potential well with sharp walls, forcing the 
particles to remain within a circle. An example of such simulation with 16 particles is shown 
in Figure [3j right. 

The first step towards considering this system as a statistical physics model is to in- 
vestigate the distribution of linear and angular velocities. For an ideal gas in 3D, the 
Maxwell-Boltzmann distribution for velocities is 




and similar for the rotational degrees of freedom. The implicit assumption in ([6]) is that the 

distribution of velocities in each direction is normal with the same width that defines the 

temperature T of the system. However, in the system of rolling particles there is no reason 

for ^ to work: first, the rotational and translational components are coupled through the 

rolling constraint and second, even though the rolling motion is along the plane, each ball 

undergoes three-dimensional motion. 

Figure [| shows several examples of the distribution for several values of total energy of 

the system. On the left side of this Figure, we plot the distribution of uj X) which is identical 

to the distribution of uj y . These distributions are always very close to normal, which are 

shown with solid curves. The right side of this Figure shows distributions of uj z for the 

same values of the energy, and it is apparent that this distribution is not normal. The 

distributions of v xi v y and v z show the same tendencies. It is important to note, however, 

that the variances (in proper units) of normally distributed angular and linear velocities 

are not the same, and thus there is no straightforward definition of temperature for this 

system. We shall also note that due to the non-normal nature of the distribution in one of 
the components and the non-holonomic coupling between angular and linear velocities, the 

kinetic energy distribution does not follow Maxwell-Boltzmann law. In addition, the kinetic 
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FIG. 4: Distributions of uj x (left) and uj z (right) for different energies of the system. Solid curves 
represent fits to normal distributions. 

energies of rotational and translational motion are not equal (even in a statistical sense), so 
the equipartition of energy in our system does not hold. 

Temperature as scaled variance Clearly, the rolling constraint prevents a straightfor- 
ward statistical physics description of the rolling particle systems. Nevertheless, there is a 
surprising relation that connects the variances a obtained from the angular (a J) and linear 
(a v ) velocities. Namely, we observed that for all values of in our numerical experiments, 
there is a surprising well-behaved linear relationship = ka V) with the coefficient k ~ 1.08 
depending only on the parameters on the system (geometry of the ball, center of mass posi- 
tion etc) but not on anything else. That relationship is valid for both lattice and gas states, 
as shown in Figure [5j We note that there is no a priori reason for such relationship to exist, 
but the surprising robustness of this law leads us to believe that it could be taken as one of 
the postulates in future development of statistical mechanics for non-holonomic gas. 




FIG. 5: Linear relationship between variances a v and for lattice states (circles) and gas states 
(lines). Common linear fit to both sets of data is also shown. The linear plots for each regime are 
shown as inserts. 

Using this surprising relation we can define the temperature as the variance of either 
horizontal linear or angular velocity, or any linear combination of those. Since the tem- 
perature is defined up to a constant, all these definitions lead to the same results up to a 
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scaling factor. Thus, on Figure [BJ we plot the scaled variances vs total energy of the system. 
If we think of the variance W as being proportional to temperature, we observe that the 
total energy E is directly proportional to the temperature, which is reminiscent of ideal gas. 
Thus, in spite of complexities of the rolling systems, some aspects from the ideal gas remain. 
We shall note that since in simulations we use real balls interacting with LJ potential, a 
more detailed simulation should show the effects of the finite size of the particles, but our 
simulations do not allow a reliable investigation of these effects. 

Using the concept of temperature as scaled invariance, we can now define thermodynam- 
ically meaningful equations of state. On Figure [BJ we plot the temperature as a function of 
the energy per particle. The equation of state for the lattice state (left panel of the Figure) 
ceases to exist for small negative energies as the lattice becomes unstable and ceases to ex- 
ist. The exact nature of the destruction of lattice state is as yet unclear, but the divergence 
of the variance a 2 ~ (E — E)~ l (shown in the insert) indicates the presence of a phase 
transition. 




E(eV) E(eV) 

FIG. 6: Variances vs energy in the lattice state (left) and gas state (right). Circles: a 2 (uj)- 
distribution; crosses: a 2 (^-distribution multiplied by k = 1.08. Left insert: l/cr 2 vs E. 

Acknowledgements The authors were partially supported by grants NSF-DMS-0908755 
and HDTRA1-10-1-0070. We also benefitted from fruitful discussions with Profs. D. D. 
Holm, P. Vorobieff and C. Tronci. 



[1] A. Ulman, Chem. Rev. 96, 1533 (1996). 

[2] S. Mamyama, T. Kimura, and M.-C. Lu, Thermal Sci. Eng. 10, 23 (2002). 
[3] R. Zangi and A. E. Mark, Phys. Rev. Lett. 91, 025502 (2003). 

[4] J. Carrasco, A. Michaelides, M. Forster, S. Haq, R. Raval, and A. Hodgson, Nat. Mater. 8, 
427 (2009). 

[5] P. Kumar, eprint arXiv:1005.2380 (2010). 

9 



[6] P. Kumar, Z. Yan, L. Xu, M. G. Mazza, S. V. Buldyrev, S. H. Chen, S. Sastry, and H. E. 

Stanley, Phys Rev Lett 97, 177802 (2006). 
[7] P. Miranda, L. Xu, Y. R. Shen, and M. Salmeron, Phys. Rev. Lett. 81, 5876 (1998). 
[8] G. Reiter, A. Demirel, and S. Granick, Science 263, 1741 (1994). 

[9] A. M. Bloch, J. E. Marsden, and D. V. Zenkov, Notices Amer. Math. Soc 52, 324 (2005). 
[10] Y. Shirai, A. J. Osgood, Y. Zhao, Y. Yao, L. Saudan, H. Yang, C. Yu-Hung, L. B. Alemany, 

T. Sasaki, J.-F. Morin, et al., J. Am. Chem. Soc. 128, 4854 (2006). 
[11] A. M. Bloch, Nonholonomic Mechanics and Control ( Springer- Verlag, 2003). 
[12] H. Goldstein, Classical Mechanics (Addison Wesley, 2001), 3rd ed. 
[13] V. E. Tarasov, Int. J. Mod. Phys. B 19, 879 (2005). 
[14] L. C. Gomes and R. Lobo, BrasRev 9, 797 (1979). 
[15] R. Kutteh, Computer Physics Communications 119, 159 (1999). 

[16] L. D. Landau and E. M. Lifshitz, Statistical Physics, part 1 (Pergamon Press, Oxford, NY, 
1980), 3rd ed. 

[17] A. Borisov and I. Mamaev, Reg. Chaot. Dyn. 7, 43 (2002). 
[18] C. G. Gray and B. G. Nickel, Am. J. Phys. 68, 821 (2000). 

[19] H. Berendsen, J. Postma, W. van Gunsteren, and J. H. in B. Pullman, ed., Reidel. Dordrecht, 

Holland p. 331 (1981). 
[20] N. M. Bou-Rabee, J. E. Marsden, and L. A. Romero, SIAM Review 50, 325 (2008). 



10 



